function res=logitRegress(X, Y, groups, nFolds)


%%%%%
% This function performs logistic regression.
%
% Inputs:
% X - input data matrix, rows are subjects and columns are covariates.
% Y - binary outcome variable, each element matching the corresponding row in X.
% nFolds - number of folds in cross-validation
%
% Outputs:
% res - structure containing logistic regression results.
%
% Written by Joon Lee, 2010
% Modified on Aug 26, 2012 by Marzyeh Ghassemi to include group balancing. 
% Modified on March 14, 2013 by Marzyeh Ghassemi to include other
%   performance specs (fscore, acc) in x-vlaidation results. 
%%%%%


res.N=size(X,1);
res.N_true=sum(Y==1);
res.N_false=sum(Y==0);

[res.b,res.dev,res.stats]=glmfit(X,Y,'binomial','link','logit');
res.yfit=glmval(res.b,X,'logit');
res.ytarget=Y;
[tpr fpr]=roc(res.ytarget',res.yfit');
res.auc=computeAUC(tpr,fpr);
res.HLtestp=HosmerLemeshowTest(res.yfit,Y);


% cross-validation    
randidx=randperm(res.N);
Xrand=X(randidx,:);
Yrand=Y(randidx);
temp=[];
yfit_xval=[];
ytarget_xval=[];
beta = [];

% cross-validation    
g0indMap = find(groups == 0);
g1indMap = find(groups == 1);
group0indices = crossvalind('Kfold', Y(groups == 0), nFolds); group0indices = group0indices(:);
group1indices = crossvalind('Kfold', Y(groups == 1), nFolds); group1indices = group1indices(:);

bestAUC = 0;

for k = 1:nFolds
    testIdx = sort([g0indMap(group0indices == k); g1indMap(group1indices == k)]);
    trainIdx= setdiff(1:length(groups), testIdx);

    b=glmfit(Xrand(trainIdx,:),Yrand(trainIdx),'binomial','link','logit');
    prob=glmval(b,Xrand(testIdx,:),'logit');
    
    yfit_xval=[yfit_xval; prob];
    ytarget_xval=[ytarget_xval; Yrand(testIdx)];
    
    % Save performance numbers
    [tpr,fpr] = roc(Yrand(testIdx)',prob');
    [sens, spec, ppv, npv, acc, fscore, str] = summaryOfPerf(Yrand(testIdx)', prob', '', 0, 0);
    temp = [temp [computeAUC(tpr,fpr); sens; spec; ppv; npv; acc; fscore]];    
    
    %Save the best beta values
    beta = [beta b];
end

res.xvalvec=temp;
res.xval=[mean(temp) std(temp)];
res.yfit_xval=yfit_xval;
res.ytarget_xval=ytarget_xval;
res.bestB = beta(:, find(temp(1, :) == max(temp(1, :))));
